function F = robustF(X, y, b, w, c)
% robustF  - compute F scores given robust stats
%
% FORMAT:       F = robustF(X, y, b, w, c)
%
% Input fields:
%
%       X           TxP model
%       y           ...xT data (n-dim supported)
%       b           ...xP beta estimates
%       w           ...xT weights (e.g. from fitrobustbisquare_img)
%       c           contrast vectors CxP for C contrasts with P weights
%
% Output fields:
%
%       F           ...xC F-contrasts
%
% Note: all values in c that are non-zero will be examined as set!

% Version:  v0.7g
% Build:    9040220
% Date:     Apr-02 2009, 8:15 PM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 5 || ...
   ~isa(X, 'double') || ...
    ndims(X) ~= 2 || ...
    size(X, 1) <= size(X, 2) || ...
    any(isinf(X(:))) || ...
    any(isnan(X(:))) || ...
   ~isa(y, 'double') || ...
   (numel(y) ~= size(X, 1) && ...
    size(y, ndims(y)) ~= size(X, 1)) || ...
    any(isinf(y(:)) | isnan(y(:))) || ...
   ~isa(b, 'double') || ...
    ndims(b) ~= ndims(y) || ...
   (numel(b) ~= size(X, 2) && ...
    size(b, ndims(b)) ~= size(X, 2)) || ...
    any(isinf(b(:)) | isnan(b(:))) || ...
   ~isa(w, 'double') || ...
    ndims(w) ~= ndims(y) || ...
    any(size(w) ~= size(y)) || ...
    any(isinf(w(:)) | isnan(w(:)) | w(:) < 0 | w(:) > 1) || ...
   ~isa(c, 'double') || ...
    ndims(c) ~= 2 || ...
    size(c, 2) ~= size(X, 2) || ...
    any(isnan(c(:)))
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing argument.' ...
    );
end

% treat single-stats case
if ndims(y) == 2 && ...
    ndims(w) == 2 && ...
    size(y, 2) == 1 && ...
    size(w, 2) == 1
    y = y';
    w = w';
end

% reshape arguments
[n, p] = size(X);
nc = size(c, 1);
sy = size(y);
sy(end) = [];
ny = prod(sy);
y = reshape(y, ny, n);
b = reshape(b, ny, p);
w = reshape(w, ny, n);
op = ones(1, p);

% simply get non-zero elements in c
ci = cell(1, nc);
for cc = 1:nc
    ci{nc} = find(c(cc, :) ~= 0);
end

% compute (raw) residual
r = y - b * X';

% compute sqrt of weights
sw = sqrt(w);

% and sum (to get idea of real d.f.)
ws = sum(w .* sw, 2);

% compute weighted residual
wr = sw .* r;
rs = sum(wr .* wr, 2);

% compute penalizing factor
pF = (1 / (n - p)) * (ws - p);

% create output data
F = zeros(ny, nc);

% test BVQXprogress
if ny > 5000
    try
        pbar = BVQXprogress;
        BVQXprogress(pbar, 'setposition', [80, 200, 640, 36]);
        BVQXprogress(pbar, 'settitle', sprintf('Computing %d contrasts...', nc));
        BVQXprogress(pbar, 0, sprintf('Running %d samples...', ny), 'visible', 0, ny);
        vcs = ceil(max(2000, ny / 100));
        vcn = vcs;
    catch
        pbar = [];
        vcn = Inf;
    end
else
    pbar = [];
    vcn = Inf;
end

% iterate over data elements (voxels)
for vc = 1:ny
    
    % create weighted matrix
    wX = sw(vc(op), :)' .* X;
    
    % iterate over contrasts
    for cc = 1:nc
        
        % compute variance explained by set elements in c
        xSS = sum((wX(:, ci{cc}) * b(vc, ci{cc})') .^ 2);
        
        % compute F test (without d.f. terms)
        F(vc, cc) = xSS / rs(vc);
    end
    
    % update progress bar
    if vc >= vcn && ...
       ~isempty(pbar)
        BVQXprogress(pbar, vc);
        vcn = vcn + vcs;
    end
end

% convert to nominal F
F(isinf(F) | isnan(F)) = 0;
for cc = 1:nc
    F(:, cc) = custom_finv(custom_fcdf( ...
        ((1 / numel(ci{cc})) * (ws - p)) .* (pF .* F(:, cc)), numel(ci{cc}), ws - p, true), ...
        numel(ci{cc}), n - p, true);
end
F = reshape(F, [sy, nc]);
F(isinf(F) | isnan(F)) = 0;

% clear progress bar
if ~isempty(pbar)
    closebar(pbar);
end
